library(specr)
library(fixest)
library(readr)
library(broom.mixed)
library(texreg)
library(dplyr)
library(fixest)
library(patchwork)

# Cross sectional linear FE models that use aggregate data (2000-2019)

source("funs_and_constants.R")

# load in data
epro_year_group_clusters <- read_csv("data/epro_data_h2.csv")

mean_h2 <- epro_year_group_clusters %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, disagreement, n_orgs, n_ed_religions, hhi_rel, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)),
    alignment_and_multiple_segments = Mode(alignment_and_multiple_segments),
    religious_constellation = Mode(religious_constellation),
    religious_favoritism_and_multiple_segments = mean(religious_favoritism_and_multiple_segments, na.rm = TRUE),
    favoritism_constellation = Mode(favoritism_constellation))


### Alignment Models


h2_bin <- feols(disagreement ~ alignment_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | countries_gwid,
                       cluster = ~countries_gwid, data = mean_h2)

# n rels as IV
h2_constellation <- feols(disagreement ~ religious_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | countries_gwid,
                              cluster = ~countries_gwid, data = mean_h2)

h2_bin_max <- feols(share_disagreement ~ alignment_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | countries_gwid,
                cluster = ~countries_gwid, data = mean_h2)

# n rels as IV
h2_constellation_max <- feols(share_disagreement ~ religious_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | countries_gwid,
                          cluster = ~countries_gwid, data = mean_h2)

screenreg(list(h2_bin, h2_constellation, h2_bin_max, h2_constellation_max), stars = stars)

##############
#### Panel Models
#############

h2_bin_fe <- feols(disagreement ~ alignment_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | year + countries_gwid,
                cluster = ~countries_gwid, data = epro_year_group_clusters)

# n rels as IV
h2_constellation_fe <- feols(disagreement ~ religious_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | year + countries_gwid,
                          cluster = ~countries_gwid, data = epro_year_group_clusters)

h2_bin_max_fe <- feols(share_disagreement ~ alignment_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | year + countries_gwid,
                    cluster = ~countries_gwid, data = epro_year_group_clusters)

# n rels as IV
h2_constellation_max_fe <- feols(share_disagreement ~ religious_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | year + countries_gwid,
                              cluster = ~countries_gwid, data = epro_year_group_clusters)


screenreg(list(h2_bin_fe, h2_constellation_fe, h2_bin_max_fe, h2_constellation_max_fe))



#################
#### logit models
#################

logit_h2 <- feglm(disagreement ~ alignment_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +  n_tek_groups | year +countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters, family = "binomial")

logit_h2_const <- feglm(disagreement ~ religious_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                              +  n_tek_groups | year + countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters, family = "binomial")


screenreg(list(logit_h2, logit_h2_const))


# make final table
texreg(list(h2_bin, h2_bin_fe, h2_constellation, h2_constellation_fe),
       stars = stars,
       custom.coef.map = coef_map,
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h2_alignment_religion_app.tex")



##################
### favoritism
################

# these are the exact same models as before but with the favoritism variable


# cross sectional
h2_bin <- feols(disagreement ~ religious_favoritism_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | countries_gwid,
                cluster = ~countries_gwid, data = mean_h2)

# n rels as IV
h2_constellation <- feols(disagreement ~ favoritism_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | countries_gwid,
                          cluster = ~countries_gwid, data = mean_h2)

h2_bin_max <- feols(share_disagreement ~ religious_favoritism_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | countries_gwid,
                    cluster = ~countries_gwid, data = mean_h2)

# n rels as IV
h2_constellation_max <- feols(share_disagreement ~ favoritism_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | countries_gwid,
                              cluster = ~countries_gwid, data = mean_h2)


screenreg(list(h2_bin, h2_constellation, h2_bin_max, h2_constellation_max))

##############
#### Panel Models
#############

h2_bin_fe <- feols(disagreement ~ religious_favoritism_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | year + countries_gwid,
                cluster = ~countries_gwid, data = epro_year_group_clusters)

# n rels as IV
h2_constellation_fe <- feols(disagreement ~ favoritism_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | year + countries_gwid,
                          cluster = ~countries_gwid, data = epro_year_group_clusters)

h2_bin_max_fe <- feols(share_disagreement ~ religious_favoritism_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | year + countries_gwid,
                    cluster = ~countries_gwid, data = epro_year_group_clusters)

# n rels as IV
h2_constellation_max_fe <- feols(share_disagreement ~ favoritism_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + n_tek_groups | year + countries_gwid,
                              cluster = ~countries_gwid, data = epro_year_group_clusters)


screenreg(list(h2_bin, h2_constellation, h2_bin_max, h2_constellation_max))



#################
#### logit models
#################

logit_h2 <- feglm(disagreement ~ religious_favoritism_and_multiple_segments + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                  +  n_tek_groups | year +countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters, family = "binomial")

logit_h2_const <- feglm(disagreement ~ favoritism_constellation + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                        +  n_tek_groups | year + countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters, family = "binomial")


screenreg(list(logit_h2, logit_h2_const))


# same table as before
texreg(list(h2_bin, h2_bin_fe, h2_constellation, h2_constellation_fe),
       stars = stars,
       custom.coef.map = coef_map,
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h2_favoritism_religion_app.tex")







